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The ability of several density-functional theory (DFT) exchange-correlation functionals to de- 
scribe hydrogen bonds in small water clusters (dimer to pentamer) in their global minimum energy 
structures is evaluated with reference to second order M0ller Plesset perturbation theory (MP2). 
Errors from basis set incompleteness have been minimized in both the MP2 reference data and the 
DFT calculations, thus enabling a consistent systematic evaluation of the true performance of the 
tested functionals. Among all the functionals considered, the hybrid X3LYP and PBEO functionals 
offer the best performance and among the non-hybrid GGA functionals mPWLYP and PBEIW 
perform the best. The popular BLYP and B3LYP functionals consistently underbind and PBE and 
PW91 display rather variable performance with cluster size. 



I. Introduction 

Density-functional theory (DFT) is the most popu- 
lar theoretical approach for determining the electronic 
structures of polyatomic systems. It has been extensively 
and successfully used to tackle all sorts of problems in 
materials science, condensed matter physics, molecular 
biology, and countless other areas. Many of these studies 
have involved the treatment of systems containing 
hydrogen bonds. Hydrogen bonds are weak (10-30 
kJ/mol ~ 100-300 meV/H bond) bonds of immense 
widespread importance, being the intermolecular force 
responsible for holding water molecules together in the 
condensed phase, the two strands of DNA in the double 
helix, and the three dimensional structure of proteins 
[l[ . A particularly important class of H-bonded systems 
are small water clusters. Small water clusters have 
been implicated in a wide range of phenomena (for 
example, environmental chemistry and ice nucleation 
di Hi) and, moreover, are thought to provide a clue as 
to the properties of liquid water. However the ability of 
DFT to quantitatively describe H bonds between H2O 
molecules in either small water clusters or the liquid 
state remains unclear. This is particularly true in light 
of recent experimental and theoretical studies which 
have raised concerns over the ability of DFT to reliably 
describe the structure and properties of liquid water 

It is now well established that the simplest approxi- 
mation to the electron exchange and correlation (XC) 
potential, the local-density approximation (LDA), is 
inappropriate for treating H bonds. For example, the 
dissociation energies of small water clusters and the 
cohesive ener gy o f ice are overestimated by >50% with 
the LDA pll. Il2l. [Tsl [T^. However, despite widespread 
practical application and several recent benchmark 
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studies it remains unclear precisely how well the many 
popular post-LDA functionals perform at describing H 
bonds between water clusters. Generalized gradient ap- 
proximation (GGA) functionals such as PBE [ll], PW91 
|la | , and BLYP flgl . for exainple, are widely used to 
examine liquid waterjlJE i, 0, ice [H [H 0, IHI 
and adsorbed water [1^, [22] , yet ask three experts which 
one is "best" and one is likely to receive three different 
answers. Likewise unanimity has not been reached 
on the performance of the many meta-GGA or hybrid 
functionals that are available, such as TPSS [11], PBEO 
HI, and B3LYP [H, [H [H, [13 . Part of the reason for 
the lack of clarity, we believe, stems from the fact that 
in previous benchmark studies insufficiently complete 
basis sets were employed and that comparisons were 
restricted to the simplest H-bonded systems involving 
H2O, namely the H2O dimer and trimer. Basis set 
incompleteness effects can, of course, mask the true 
performance of a given functional and, as we will show 
below, the ability of a given functional to accurately 
predict the strength of the H bond in the dimer or even 
the trimer does not necessarily reveal how well that 
functional will perform even for the next largest clusters, 
tetramers and pentamers. 



In the following we report a study in which the ability 
of several GGA, meta-GGA, and hybrid functionals to 
compute the energy and structure of H bonds between 
H2O molecules is evaluated. So as to enable the use 
of large basis sets, which we demonstrate approach the 
complete basis set (CBS) limit, in the generation of the 
benchmark data and the DFT data itself, this study is 
limited to the four smallest H2O clusters (dimer, trimer, 
tetramer, and pentamer). In addition, this study is 
restricted to the established lowest energy conformer of 
each cluster [ll,[il, [13], which, for orientation purposes, 
we show in Fig. 1. For this admittedly small structural 
data set we find that, of the functionals tested, the 
hybrid X3LYP [U and PBEO [^ functionals offer the 
best performance. Among the regular (pure) GGAs 
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mPWLYP [H,!!! and PBEIW perform best. BLYP 
[H [11 and B3LYP [H, [H HE |27i predict too weak H 
bonds and PEE and PW91 |l6l | display rather vari- 
able performance with cluster size. Although MPWBIK 
[Hi, PW6B95 35], and B98 [H] were previously shown 
to offer outstanding performance for water, this is not 
now the case when highly accurate basis sets are used. 

II. Reference data - MP2 

For a systematic benchmark study such as this, 
reliable reference data is essential. Experiment is, in 
principle, one source of this data. However, experi- 
mental dissociation energies are simply not available 
or do not come with sufhciently small error bars for 
all the H2O clusters we examine here. Further, with 
our aim to systematically evaluate the performance of 
many DFT XC functionals it becomes impractical to 
compute all the small contributions to the experimental 
dissociation energy that come on top of the total 
electronic dissociation energy - an easily accessible 
total energy difference - such as zero point vibrations, 
relativistic contributions, etc. The obvious alternative 
source of reference data are the results obtained from 
correlated quantum chemistry methods such as second 
order M0ller Plesset perturbation theory (MP2) [s^l 
or coupled-cluster theory [sl]. Indeed such methods 
have been widely a.pplied to examine H-bonded systems 

[n, m, M mm\m, \m, m\m,m\m m m m- 

In particular coupled-cluster with single and double 
excitations plus a perturbative correction for connected 
triples (CCSD(T)) produces essentially "exact" answers 
if sufhciently accurate basis sets are used. For example, 
the best CCSD(T) value for the binding energy of the 
water dimer is at 217. 6±2 meV [51] in good agreement 
with the appropriate experimental number of 216.8±30 
meV [52, J53j. However, since the computational cost 
of CCSD(T) formally scales as N'^ , where N is the 
number of basis functions, the most extravagant use of 
computational power is required for CCSD(T) calcula- 
tions with large basis sets. MP2, on the other hand, 
scales as and when compared to CCSD(T) for water 
dimers and trimers at the CBS limit, yields binding 
energies that differ by no more than 2 meV/H bond 
[43I . 0| ■ In addition, a recent study of water hexamers 
using CCSD(T) with an aug-cc-pVTZ basis set revealed 
that the MP2 and CCSD(T) dissociation energies of 
various hexamer structures differ by <3 mcV/H20 [s^ . 
Thus MP2 is a suitable method for obtaining reference 
data with an accuracy to within a few meV/H bond. 
Such accuracy, which is well beyond so-called chemical 
accuracy (Ikcal « 43 meV), is essential in studies of 
H-bonded systems. 

Since MP2 geometries are not available for all four 
clusters examined here we have computed new MP2 
structures for each one. All calculations have been 
performed with the Gaussian03 [s^ and NWChem 
[56| codes and all geometries were optimized with 
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FIG. 1: Structures of the four water clusters examined here 
in their global minimum energy configurations. The dashed 
lines indicate H bonds. Some of the structural parameters of 
the H bond are indicated alongside the dimer. We note that 
in the trimer, tetramer, and pentamer there is one H bond 
per water molecule. 



an aug-cc-pVTZ basis set within the "frozen core" 
approximation i.e., correlations of the oxygen Is orbital 
were not considered [H^- Although the aug-cc-pVTZ 
basis set is moderately large (92 basis functions/H20), 
this finite basis set will introduce errors in our predicted 
MP2 structures. However, a test with the H2O dimer 
reveals that the aug-cc-pVTZ and aug-cc-pVQZ MP2 
structures differ by only 0.004 A in the 0-0 bond length 
and 0.16° in the H bond angle {(j>, Fig. 1). Likewise, 
Nielsen and co-workers have shown that the MP2 0-0 
distances in the cyclic trimer differ by 0.006 A between 
the aug-cc-pVTZ and aug-cc-pVQZ basis sets with all 
other bonds differing by <0.003 A ^49^. For our present 
purposes these basis set incompleteness errors on the 
structures are acceptable and it seems reasonable to 
assume that the aug-cc-pVTZ structures reported here 
come with error bars of ±0.01 A for bond lengths and 
±0.5° for bond angles. 

Total energies and dissociation energies are known to 
be more sensitive to basis set incompleteness effects than 
the geometries are. To obtain reliable MP2 total en- 
ergies and dissociation energies we employ the aug-cc- 
pVTZ, aug-cc-pVQZ (172 basis functions/H20) and aug- 
cc-pV5Z (287 basis functions/H20) basis sets in conjunc- 
tion with the well-established methods for extrapolating 
to the CBS limit. Usually the extrapolation schemes rely 
on extrapolating separately the Hartree-Fock (HF) and 
correlation contributions to the MP2 total energy. For 
extrapolation of the HF part we use Feller's exponential 
fit [11]: 
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FIG. 2: (a) Variation in the MP2 dissociation energy for tlie 
H2O dimer without a counterpoise correction for basis set 
superposition error (BSSE) (labeled MP2) and with a coun- 
terpoise correction for BSSE (labeled MP2(CP)) as a function 
of basis set size. The extrapolated complete basis set (CBS) 
dissociation energy for the H2O dimer with MP2 is also indi- 
cated, (b) Variation in the dissociation energy for the H2O 
dimer with and without a counterpoise BSSE correction as a 
function of basis set size for three different DFT functionals. 
The basis set labels on the X axis of (a) and (b) indicate aug- 
cc-pVXZ basis sets, where X=3, 4, and 5. Lines are drawn 
to guide the eye only. All structures were optimized with an 
aug-cc-pVTZ basis set consistently with MP2 and with each 
DFT functional. 



basis set {X—3, 4, and 5 for the aug-cc-pVTZ, aug-cc- 
pVQZ, and aug-cc-pV5Z basis sets, respectively). Ei^^ is 
the corresponding HF energy, E^gg is the extrapolated 
HF energy at the CBS limit, and A and B are fitting 
parameters. For the correlation part of the MP2 total 
energy we follow an inverse power of highest angular mo- 
mentum equation [s^ [g^, [61| ■ 
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where E^°^^ is the correlation energy corresponding to 
X, E^'g'g is the extrapolated CBS correlation energy, 
and C and D are fitting parameters. We have tested 
various extrapolation schemes available in the literature 
[58|, I59i led, l6l|, ISi l63j, Q and did not see a differ- 
ence of more than 1.2 meV/H bond between all the pre- 
dicted CBS values. We opted for the scheme provided by 
eqns. ([T][2]) because we found that with input from triple- 
, quadruple-, and pentuple-C basis sets this method was 
best able to predict the total energy of a water monomer 
and dimer explicitly calculated with an aug-cc-pV6Z ba- 
sis set (443 basis functions/H20). Having obtained MP2 
CBS total energies for the H2O monomer and each of the 
H2O clusters, we thus arrive at the MP2 CBS electronic 
dissociation energies (f") per H bond which are given 
by, 
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where E^^^^ is the total energy of each cluster with 
n H2O molecules, E^^'^ is the total energy of a H2O 
monomer, and ?^//— bond 

is the number of H bonds in 



the cluster. Our CBS MP2 binding energies for the 
dimer, trimer, tetramer, and pentamer are 215.8, 228.5, 
299.9, and 314.4 meV/H bond, respectively These 
values are all within 0.5 meV/H bond of the previous 
MP2 CBS dissociation energies reported by Xantheas et 
al. [iOI. We expect that the various errors accepted in 
producing these values (MP2 (valence only) treatment 
of correlation, aug-cc-pVTZ structures, extrapolation to 
reach the CBS, etc.) will lead to errors in our reference 
data from the exact electronic dissociation energies on 
the order of ±5.0 meV/H bond at most. With our 
present aim to evaluate the performance of various DFT 
XC functionals such errors are acceptable. 

III. DFT 

In a study such as this there is an essentially endless 
list of functionals that we could consider evaluating. 
Here we have chosen to examine 16 different func- 
tionals which are widely used or have been reported 
to perform particularly well for H-bonded systems in 
predicting dissociation energies and structures of the 
above mentioned clusters. Specifically we have chosen 
to optimize structures of each cluster with the following 
functionals: (I) PW91 [l^ - an extremely popular 
non-em piri cal GGA widely used in calculations of bulk 
ice [Iljel, and other H-bonded systems [3; (II) 
PBE [is'l - the twin of PW91 that has again been widely 
used and tested for H-bonded systems [H, ; (III) 

PBEIW - a parameterized empirical variant of PBE 
specifically designed to yield improved energetics of H 
bonds [33]. (IV) TPSS 23] - the meta-GGA variant of 
PBE, recently used in simulations of liquid water and 
evaluated for small water clusters [I,!!!!!!!!!!!]; (V) 
PBEO 241] - a so-called parameter free hybrid variant 
of PBE, also recently tested for water 0, [H, IH, ill; 
(VI) BLYP - Becke88 jT^I exchange combined with LYP 
|18| correlation, a popular functional for liquid water 
simulations i,i,i,0,[ll; (VII) B3LYP [H, [M S [13 
- the extremely popular Beckc three parameter hybrid 
functional combined with LYP nonlocal correlation, 
which has, of course, been widely used to examine 
H-bonded systems |4, 42, il S^; (VIII) mPWLYP - 
a combination of a modified PW91 exchange functional 
(mPW) [13 with the LYP correlation functional, found 
to be the most accurate pure GGA for the energetics of 
H bonds in water dimers and trimers [33] . (IX) BP86 - 
an empirical GGA combining Becke88 y/Tj] exchange and 
Perdew86 [t^ correlation that is well-tested for hydrogen 
bonded systems [is^, (X) X3LYP [H - another 

empirical hybrid functional designed to describe weak 
(non-covalent) interactions that is becoming a familiar 
name for calculations of water [1, [H, [Hol; (XI) XLYP 
[33 - the non-hybrid GGA version of X3LYP, also tested 
for H-bonded systems [13; (XII) B98 0] - another 
hybrid functional, said to perform extremely well for 
water clusters [11, [13; (XIII) MPWBIK Jp - a one 
parameter hybrid meta-GGA using mPW [32] exchange 
and Becke95 [73 correlation, said to be the joint-best 
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for H bonds between water molecules [H, (XIV) 
PW6B95 [111 - another hybrid meta-GGA combining 
PW91 [IB] exchange and Becke95 [lH correlation, found 
to be the other joint-best functional for the H bonds 
between water molecules (33j; (XV) B3P86 - Becke 3 
parameter hybrid functional combined with Perdew86 
nonlocal correlation, found to be best functional 
for H-bonded systems in a recent benchmark study 
[12; and (XVI) BH&HLYP [13, [H, [zl - said to offer 
similar performance to B3P86 for H-bonded systems [i^ . 

As with MP2, the question arises as to what basis 
sets to use in order to ensure that the DFT results 
reported here are not subject to significant basis set 
incompleteness errors, which would cloud our evaluations 
of the various functionals. There are no established ex- 
trapolation schemes for DFT. However, it is well-known 
that DFT total energies are less sensitive to basis set 
size than explicitly correlated methods such as MP2 
[li [li, lli, [li- Indeed from the plot in Fig. 2 it can 
be seen that the computed DFT dissociation energies 
converge much more rapidly with respect to basis set size 
than MP2 does (c./. Figs. 2(a) and (b)). Specifically we 
find that upon going from aug-cc-pVTZ to aug-cc-pV5Z 
the dissociation energy of the H2O dimer changes by only 
1.0, 2.7, and 1.5 meV for the PBE, TPSS, and PBEO 
functionals, respectively. Further, with the aug-cc-pV5Z 
basis set we find that the counterpoise corrected and 
uncorrected dissociation energies essentially fall on top 
of each other, with the largest difference for the dimers 
and trimers being 0.45 meV/H bond with the TPSS 
functional. In addition, upon going beyond aug-cc-pV5Z 
to aug-cc-pV6Z the dimcr dissociation energies change 
by only 0.24, 0.11, 0.19, 0.25 meV for the PBE, TPSS, 
PBEO and BLYP functionals, respectively. Thus the 
DFT dissociation energies we report in the following 
will all come from those obtained with the aug-cc- 
pV5Z basis set, which is sufficiently large to reflect 
the true performance of each functional at a level of 
accuracy that is reasonably expected to approach the 
basis set limit to within about 0.5 meV/H bond or better. 

IV. Results 

A. Dissociation energy 

In Table [T] the computed dissociation energies obtained 
with MP2 and with each of the DFT functionals are re- 
ported. To allow for a more convenient comparison of 
the performance of the various functionals we plot in Fig. 
ini^a) the difference between the DFT and MP2 dissocia- 
tion energies (AD") as a function of water cluster size. In 
this figure positive values correspond to an overestimate 
of the dissociation energy by a given DFT functional com- 
pared to MP2. So, what do we learn from Table 1 and 
Fig. [21Ia)? First, the functionals which offer the best per- 
formance for the clusters examined are the hybrid X3LYP 
and PBEO functionals, coming within 7 meV/H bond 
for all four clusters. Of the non-hybrid functionals the 
pure GGAs mPWLYP and PBEIW perform best, com- 



ing within 12 meV/H bond for all four clusters. Second, 
the very popular BLYP and B3LYP functionals consis- 
tently underbind: B3LYP predicts H bonds that are '^20 
meV too weak; and BLYP predicts H bonds that are '^35 
meV too weak. Third, PBE overestimates the binding in 
the dimer and trimer ever so slightly, coming within 5 
meV/H bond, but for the tetramer and pentamer drifts 
away to yield errors of ^--^20 meV/H bond. Fourth, PBE 
and PW91 exhibit a non-negligible difference. Although 
it is often assumed that identical numerical results should 
be obtained from these two functionals this is not the 
case here; PW91 is consistently 12-14 meV/H bond worse 
than PBE. Both functionals, however, exhibit a similar 
tendency towards increased overbinding as the cluster 
size grows. Indeed it is clear from Fig. UJa) that all 
PBE-related functionals (PBE, PW91, PBEIW, TPSS, 
and PBEO) show this trend, which in the case of TPSS 
means that it gets within ^7 meV/H bond for the pen- 
tamer starting from an error of '^20 meV/H bond for 
the dimer. Likewise PBEIW gets closer to the reference 
value as the cluster size grows. Finally, despite previous 
suggestions to the contrary [3l,[42; 69], none of the other 
functionals particularly stand out: B98 underbinds by 
just over 13 meV/H bond, and BP86 exhibits a rather 
strong variation in performance with cluster size, rang- 
ing from a 30 to 14 mcV/H bond error. B3P86 shows 
similar behavior to BP86, although the magnitude of the 
error is much less and indeed for the tetramer and pen- 
tamer B3P86 gives values close (within 3 meV/H bond) 
to MP2. MPWBIK and PW6B95 both underbind by 
>20 meV/H bond. 

B. Cooperativity 

An important aspect of the energetics of H bonds is 
that they tend to undergo cooperative enhancements, 
which for the present systems implies that the average 
strengths of the H bonds between the water molecules 
increases as the number of H bonds increases. The fact 
that the H bonds in water clusters undergo cooperative 
enhancements is now well established [H, [2 [Z^] , as is the 
importance of cooperativity in many other types of H- 
bonded systems [l[ [12, [t^] . Here we have evaluated the 
ability of each functional to correctly capture the com- 
puted MP2 cooperative enhancement, defined as the per- 
centage increase in the average H bond strength com- 
pared to that in the H2O dimer. These numbers are 
reported in parenthesis in Table [J We find that all func- 
tionals capture the correct trend, i.e., the average H bond 
strength increases upon going from dimer to pentamer. 
In addition, most functionals get the absolute percentage 
enhancement correct to within 5%. The notable excep- 
tions are BP86, B3P86, and TPSS which for the tetramer 
and pentamer predict cooperative enhancements that ex- 
ceed the MP2 values by 10-15%. 

C. Geometry 

Let us turn now to an assessment of the quality of the 
structural predictions made by each DFT functional. 
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TABLE I: Comparison of the MP2 complete basis set dissociation energies to those obtained with various DFT functionals 
computed with an aug-cc-pV5Z basis set for four different water clusters. DFT dissociation energies that come within ±5.0 
meV of the corresponding MP2 value are indicated in bold. The numbers in parenthesis indicate the percentage cooperative 
enhancement in the H bond strength compared to the dissociation energy of the dimer. Averages of the signed and unsigned 
errors in the dissociation energies of all DFT functionals from the corresponding MP2 numbers over all four clusters are also 
provided as ME (moan error) and MAE (mean absolute error). The DFT functionals are ordered in terms of increasing MAE. 
All structures wore optimized consistently with MP2 and with each DFT functional with an aug-cc-pVTZ basis set and all 
values are in meV/H bond (IKcal/mol = 43.3641 meV). 





Dimer 


Trimer 


Tetramer 


Pentamer 


ME 


MAE 


MP2 


215.8 


228.5 (5.9) 


299.9 (38.9) 


314.4 (45.7) 






X3LYP 


213.8 


221.9 (3.8) 


298.3 (39.5) 


316.0 (47.8) 


-2.2 


2.9 


PBEO 


214.5 


224.6 (4.7) 


302.7 (41.1) 


320.9 (49.6) 


1.0 


3.6 


mPWLYP 


218.5 


226.0 (3.4) 


305.4 (39.8) 


323.7 (48.1) 


3.8 


5.0 


B3P86 


203.5 


220.0 (8.1) 


299.4 (47.1) 


316.5 (55.5) 


-4.8 


5.9 


PBEIW 


207.9 


216.6 (4.0) 


294.9 (41.8) 


312.7 (50.4) 


-6.6 


6.6 


BH&HLYP 


213.2 


219.5 (2.9) 


291.3 (36.6) 


308.3 (44.6) 


-6.6 


6.6 


PBE 


220.1 


233.5 (6.1) 


316.4 (43.8) 


334.8 (52.1) 


11.6 


11.6 


B98 


205.6 


211.4 (2.8) 


285.9 (39.1) 


303.1 (47.4) 


-13.2 


13.2 


TPSS 


196.4 


209.4 (6.6) 


288.8 (47.0) 


307.5 (56.6) 


-14.1 


14.1 


B3LYP 


197.4 


206.3 (4.5) 


280.1 (41.9) 


297.2 (50.6) 


-19.4 


19.4 


PW6B95 


200.9 


210.5 (4.8) 


276.8 (37.8) 


292.7 (45.7) 


-19.4 


19.4 


MPWBIK 


199.1 


210.6 (5.5) 


276.3 (38.8) 


292.3 (46.8) 


-20.1 


20.1 


BP86 


184.4 


205.7 (11.6) 


282.5 (53.2) 


300.8 (63.1) 


-21.3 


21.3 


PW91 


232.5 


244.9 (5.1) 


330.8 (42.3) 


350.5 (50.8) 


25.0 


25.0 


XLYP 


191.4 


198.6 (3.8) 


272.2 (42.2) 


288.9 (50.9) 


-26.9 


26.9 


BLYP 


180.7 


191.7 (6.1) 


264.9 (46.6) 


281.2 (55.6) 


-35.0 


35.0 




FIG. 3: (a) Difference in the dissociation energy (AD") in meV/H bond of the various DFT functionals compared to MP2, 
plotted as a function of cluster size. Positive values correspond to an overestimate of the dissociation energy by a given DFT 
functional, (b) Average value of the MP2 and DFT 0-0 distances (Ro-o) as a function of cluster size. The inset zooms in 
on the dimer region, (c) Difference in the average O-O distance (ARo-o) between MP2 and DFT. Positive values correspond 
to an overestimate of the average 0-0 distances by a given DFT functional, (a)-(c) All DFT energies are calculated with an 
aug-cc-pV5Z basis set on geometries optimized consistently with each functional with an aug-cc-pVTZ basis set. Lines are 
drawn to guide the eye only. 
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TABLE II: Mean absolute error (MAE) of various DFT functionals from MP2 for five different structural parameters, averaged 
over the four water clusters examined here. The numbers in bold all have MAE <0.010 A for bond lengths and <0.50° for 
bond angles. Mean errors (ME) are given in parenthesis. All structures were optimized consistently with MP2 and with each 
DFT functional with an aug-cc-pVTZ basis set. The order of the functionals is the same as in Table I. 





ARo-o (A) 


ARhb (A) 


ARo-H (A) 


Acp (°) 


Ae (°) 


X3LYP 


0.002 (-0.002) 


0.003 (-0.003) 


0.001 (0.000) 


0.21 (0.21) 


1.04 (1.04) 


PBEO 


0.024 (-0.024) 


0.023 (-0.023) 


0.002 (-0.001) 


0.77 (0.77) 


0.69 (0.69) 


mPWLYP 


0.012 (0.012) 


0.008 (-0.004) 


0.012 (0.012) 


0.61 (0.47) 


0.51 (0.51) 


B3P86 


0.042 (-0.042) 


0.051 (-0.051) 


0.003 (0.001) 


1.00 (1.00) 


0.77 (0.77) 


PBEIW 


0.011 (0.009) 


0.010 (-0.006) 


0.011 (0.011) 


1.13 (1.13) 


0.13 (0.13) 


BH&HLYP 


0.006 (-0.003) 


0.015 (0.015) 


0.013 (-0.013) 


0.48 (-0.17) 


1.52 (1.52) 


PBE 


0.024 (-0.024) 


0.046 (-0.046) 


0.012 (0.012) 


1.43 (1.21) 


0.13 (0.13) 


B98 


0.016 (0.016) 


0.015 (0.015) 


0.001 (-0.001) 


0.52 (0.52) 


0.66 (0.66) 


TPSS 


0.018 (-0.018) 


0.037 (-0.037) 


0.010 (0.010) 


1.28 (1.25) 


0.22 (0.22) 


B3LYP 


0.009 (0.009) 


0.007 (0.007) 


0.001 (0.001) 


0.31 (0.31) 


0.93 (0.93) 


PW6B95 


0.026 (0.026) 


0.029 (0.029) 


0.006 (-0.006) 


0.29 (0.24) 


0.81 (0.81) 


MPWBIK 


0.006 (0.006) 


0.016 (0.016) 


0.012 (-0.012) 


0.38 (0.31) 


1.09 (1.09) 


BP86 


0.028 (-0.028) 


0.051 (-0.051) 


0.014 (0.014) 


1.58 (1.46) 


0.11 (0.11) 


PW91 


0.038 (-0.038) 


0.038 (-0.038) 


0.012 (0.012) 


1.44 (1.21) 


0.29 (0.29) 


XLYP 


0.040 (0.040) 


0.028 (0.028) 


0.011 (0.011) 


0.53 (0.49) 


0.37 (0.37) 


BLYP 


0.031 (0.031) 


0.015 (0.015) 


0.009 (0.009) 


0.69 (0.64) 


0.37 (0.37) 



The five key structural parameters of the H2O clusters 
that we evaluate are: (i) The distance between adjacent 
oxygen atoms involved in a H bond, Ro-o; (ii) The 
length of a H bond, given by the distance between the 
donor H and the acceptor O, Rq - h — Rhb (Fig. 1); (iii) 
The H bond angle, Z(0 • • • H-0) = (Fig. 1); (iv) The 
internal 0-H bond lengths of each water, Rq-h; and (v) 
The internal H-O-H angle of each water, Z(H-O-H) — 9 
(Fig. 1). In Table II the mean absolute error (MAE) 
and mean error (ME) between the MP2 values and those 
obtained from each functional, averaged over all four 
clusters, are listed for each of the above parameters. This 
provides an immediate overview for how the functionals 
perform. Summarizing the results of this table, we 
find that X3LYP, BH&HLYP, B3LYP, and MPWBIK 
perform the best for 0-0 distances. All those function- 
als yield results that are essentially identical to MP2, 
coming within our estimated MP2 bond distance error 
bar of 0.01 A. B3P86 is the worst functional in terms 
of 0-0 distances, with a MAE of 0.04 A. Largely, these 
conclusions hold for the related quantity, Rhb, although 
now B3P86, BP86, and PBE perform worst with MAE 
values of ^0.05 A. In terms of the H bond angle, 0, 
X3LYP, B3LYP, PW6B95, MPWBIK, and BH&HLYP 
are essentially identical to MP 2 coming within our 
estimated MP 2 error bar for angles of 0.5° and again 
PW91, PBE, and BP86 are the worst being -1.5° 
away from MP2. For the internal 0-H bond lengths no 
functional is worse than ^^0.015 A and for the internal 
H-O-H angles, 9, all functionals are within —1.5° of MP2. 

One specific aspect of the structures of the small 
cyclic water clusters examined here, that is known from 
experiment and previous calculations 0, 1?^ is that the 
average 0-0 distances between the H2O molecules in the 



clusters shorten as the cluster size increases. This trend 
is, of course, related to the cooperative enhancement 
in H bond strengths discussed already. As can be seen 
from the plot of computed 0-0 distances versus cluster 
size in Fig. all functionals correctly capture this 

effect: the ~0.2 A shortening in the 0-0 bond distances 
upon going from dimer to pentamer predicted by MP2 
is also captured by every DFT functional. To look at 
this issue more closely and specifically to examine how 
each functional varies with respect to MP2 we plot in 
Fig. 3(c) the difference between the MP2 and DFT 0-0 
distances for the four clusters. Positive values in Fig. 
ISlJc) indicate that the DFT 0-0 bonds are longer than 
the MP2 ones. We note that the average MP2 0-0 
distances for the dimer, trimer, tetramer and pentamer 
are 2.907, 2.787, 2.732, and 2.716 A, respectively. As 
indicated already in our previous discussion, X3LYP, 
B3LYP, BH&HLYP, and MPWBIK perform the best at 
predicting the correct 0-0 bond length for each cluster; 
coming within 0.01 A of the MP2 values on every 
occasion. Indeed the consistent closeness of the X3LYP 
0-0 distances to the MP2 ones is remarkable. PBEO 
is a little worse than X3LYP for the 0-0 distances, 
predicting bonds which are consistently about 0.02-0.03 
A too short. Of the other functionals B3P86 stands out 
as predicting the shortest 0-0 distances (always —0.04 
A less than MP2) and XLYP and BLYP predict the 
longest 0-0 distances, always at least 0.02 A longer 
than MP 2. 

V. Discussion 

Here we have shown that of the functionals tested 
X3LYP and PBEO offer exceptional performance for the 
H bonds in small water clusters in their global minimum 
energy structures. However, a previous benchmark study 
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on the ability of most of the functionals considered here 
to describe the energetics of H bonds between water 
molecules has arrived at somewhat different conclusions 
dl]. Specifically, a MAE of 19.5 meV/H bond has 
been reported for PBEO, worse than the MAE of 3.6 
meV/H bond obtained here. In addition, MAEs of 5-7 
meV/H bond have been reported with the PW6B95, 
MPWBIK, and B98 functionals, suggesting improved 
performance for these functionals over what we find here. 
In that study the so-called MG3S basis set (identical 
to 6-311-|-G(2df,2p) for H2O) was used. By comparing 
the performance of the above-mentioned functionals 
with the MG3S and the aug-cc-pV5Z basis sets for 
the four clusters under consideration here it appears 
that the incompleteness of the MG3S basis set is the 
main reason for the small discrepancy. The results, 
illustrated in the histogram in Fig. IH reveal that the 
dissociation energies obtained with the MG3S basis 
set are consistently ~18 meV (0.42 kcal/mol) per H 
bond larger than those obtained with the aug-cc-pV5Z 
basis set. Thus ahhough PW6B95, MPWBIK, and B98 
perform well with the MG3S basis set (all within ±7 
meV/H bond of MP2 for the clusters considered here) 
all exhibit a propensity to underbind when the more 
complete aug-cc-pV5Z basis set is used. Conversely, 
PBEO and one other functional tested, mPWLYP, which 
predict too strong H bonds with the MG3S basis set 
(MAEs of 18.1 and 22.3 meV/H bond for the PBEO 
and mPWLYP functionals, respectively, for the clusters 
examined here) actually perform very well with the more 
complete aug-cc-pV5Z basis set (MAEs of 3.6 and 5.0 
meV/H bond for the PBEO and mPWLYP functionals, 
respectively). The small and systematic overbinding due 
to the incompleteness of the MG3S basis set has also 
been pointed out by Csonka et al. [68j . 

Another interesting aspect of the results of the present 
study is that the performance of some functionals differs 
appreciably from one cluster to another. For example, 
PBE is only ~4-5 meV/H bond away from MP2 for 
the dimer and trimer but >15 meV/H bond away 
from MP2 for the tetramer and pentamer. Conversely, 
TPSS is -20 meV/H bond off MP2 for the dimer but 
within 7 meV/H bond of MP2 for the pentamer. Other 
functionals which show strong variation in performance 
with cluster size are PW91, BP86, and B3P86, and 
the functional in the admirable position of showing the 
least variation, consistently predicting H bonds that are 
^--^35 meV too weak, is BLYP. The general conclusion 
of this analysis, however, is that it is not necessarily 
sufhcient to use the performance of a given functional 
for a single system, such as for example the H2O dimer, 
as a guide to how that functional will perform for H 
bonds between H2O molecules in general. Indeed the 
results reported here indicate that H bond test sets such 
as the "W7" test set [l^l for water would benefit from 
the inclusion of structures other than dimers and trimers. 
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FIG. 4: Mean error (ME) in the dissociation energies obtained 
with aug-cc-pV5Z (solid bars) and MG3S (hashed bars) basis 
sets for five selected functionals for the four clusters examined 
here. Positive values correspond to an average overestimate 
of the dissociation energy compared to MP2 for the clusters. 
All errors are measured relative to our reference CBS MP2 
values. 



We now ask if the results and conclusions arrived at 
here are of general relevance to H2O molecules in other 
environments and to other types of H-bonded systems. 
Some parallels with DFT simulations of liquid H2O 
can be seen. It is generally found, for example, that 
(when everything else is equivalent) BLYP liquid H2O 
is less structured (i.e., the first peak of the 0-0 radial 
distribution function (RDF) has a lower maximum) 
than PBE hquid H2O ;4, J, H, 0, 111; consistent with 
the weaker H bonds predicted by BLYP compared to 
PBE. Similarly, the first simulations of liquid H2O 
with hybrid DFT functionals (B3LYP, X3LYP, and 
PBEO) have recently been reported Q and the trend 
in the position of the first peak in the 0-0 RDF can 
be interpreted as being consistent with the current 
observations. Specifically it was found (although the 
error bars are large because the simulations were short (5 
ps)) that the position of the first peak in the 0-0 RDF 
moves to shorter separation upon going from B3LYP 
to X3LYP to PBEO, which is consistent with the small 
decrease of the 0-0 distances (Fig. WL^)) and increase 
in H bond strengths along this series (Table I). Looking 
at other H-bonded systems with slightly stronger (for 
example, NH3 • • -1120) or slightly weaker H bonds (for 
example, NH3 • • -NHs) than those considered here it is 
known, for example, that PBE generally overestimates 
these H bond strengths slightly: PBE overestimates 
NH3 • • ■ H2O by -30 meV and NH3 • • • NH3 by 6 meV 
[4^. Likewise, BLYP and B3LYP have been shown to 
underestimate a range of H-bonded systems by 20-30 
meV/H bond [46l |. However, the general performance 
of X3LYP and PBEO for other H bonded systems has 
not been evaluated yet in any great detail with suitably 
large basis sets. In light of the present results it will 
be interesting to see how well these functionals perform 
for other H-bonded systems. Likewise mPWLYP and 



8 



PBEIW are not widely used. Since they are pure GGAs 
(without any contribution from HF exchange) they will 
offer computational savings compared to X3LYP and 
PBEO, particularly for condensed phase simulations, and 
would thus be interesting to explore further for other 
H-bonded systems. 

Finally, we point out that an interesting conclusion 
of the present study is the non-negligible difference 
between the H bond energies predicted by PBE and 
PW91; with PW91 consistently being 12-14 meV/H 
bond worse than PBE. A similar discrepancy, although 
in a rather different area of application - surface and 
defect formation energies of metals - has been identified 
by Mattsson and co-workers fs^. Specifically they 
found that the PW91 and PBE monovacancy formation 
energies of Al differed by ^--^30-40 meV. As Mattson and 
co-workers have done, we caution that it does now not 
seem wise to expect identical numerical results from 
PBE and PW91. 

VI. Conclusions 

In summary, we have computed MP 2 CBS values for 
the dissociation energies of small H2O clusters (dimer to 
pentamer) in their global minimum energy structures. 
This data has been used to evaluate the performance 
of 16 DFT functionals. All DFT energies reported here 
have been obtained with an aug-cc-pV5Z basis set, 
which for DFT is sufficiently large to enable the true 
performance of each functional to be assessed, absent 
from significant basis set incompleteness errors. Among 
the functionals tested we find that PBEO and X3LYP 
perform best for the energetics of the H bonds considered 
here; always being within 10 mcV/H bond of MP2. 
In terms of the structures X3LYP offers outstanding 
performance, predicting structures essentially identical 
to MP2 for all four clusters. Of the pure GGAs con- 
sidered mPWLYP and PBEIW perform best. A small 



but non-negligible difference in the results obtained 
with PBE and PW91 has been identified, with PBE 
consistently being 12-14 meV/H bond closer to MP2 
than PW91. 

In closing we note that, although X3LYP and PBEO 
predict the most accurate H bond energies, it is impor- 
tant to remember that all functionals considered here 
do reasonably well. If, for example, one's definition of 
"good" is so-called chemical accuracy (Ikcal/mol « 43 
meV/H bond) then it is clear from Fig. [3Ja) that all 
functionals achieve chemical accuracy for all clusters. 
The problem is, of course, that for bonds as weak as H 
bonds, chemical accuracy is a rather loose criterion since 
it amounts to around 20-30% of the total bond strength. 
Future work will involve the investigation of larger 
H2O clusters in which the ability of DFT functionals to 
correctly describe the ordering of low energy isomeric 
structures becomes crucial. 
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basis set) of all the cluster studied here in xyz format. 
The total energies of each clusters obtained from MP2 
and the 16 DFT functionals are also provided. This 
document can be reached through . 
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